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ABSTRACT 

Within the area of systems health management, the task of 
prognostics centers on predicting when components will fail. 
Model-based prognostics exploits domain knowledge of the 
system, its components, and how they fail by casting the un- 
derlying physical phenomena in a physics-based model that is 
derived from first principles. Uncertainty cannot be avoided 
in prediction, therefore, algorithms are employed that help in 
managing these uncertainties. The particle filtering algorithm 
has become a popular choice for model-based prognostics due 
to its wide applicability, ease of implementation, and support 
for uncertainty management. We develop a general model- 
based prognostics methodology within a robust probabilistic 
framework using particle filters. As a case study, we consider 
a pneumatic valve from the Space Shuttle cryogenic refuel- 
ing system. We develop a detailed physics-based model of 
the pneumatic valve, and perform comprehensive simulation 
experiments to illustrate our prognostics approach and evalu- 
ate its effectiveness and robustness. The approach is demon- 
strated using historical pneumatic valve data from the refuel- 
ing system. 

1 Introduction 

Prognostics is concerned with determining the health of sys- 
tem components and making end of life (EOL) and remain- 
ing useful life (RUL) predictions. It is a key enabling tech- 
nology for condition-based maintenance, and serves to in- 
crease system availability, reliability, and safety by enabling 
timely maintenance decisions to be made. As with diagnos- 
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tics, prognostics methods are typically categorized as either 
model-based or data-driven. Data-driven approaches, rather 
than taking advantage of system and domain knowledge, in- 
stead utilize large amounts of run-to-failure data that are used 
to train machine learning algorithms to identify trends and 
determine EOL and RUL (Schwabacher, 2005). However, 
such data is often difficult to acquire. In contrast, model- 
based approaches exploit domain knowledge of the system, 
its components, and how they fail, in the form of models, in 
order to provide EOL and RUL predictions (Roemer, Bying- 
ton, Kacprzynski, & Vachtsevanos, 2005; Byington, Watson, 
Edwards, & Stoelting, 2004; Saha & Goebel, 2009; Daigle & 
Goebel, 2010b; Luo, Pattipati, Qiao, & Chigusa, 2008). The 
underlying physical phenomena are captured in a physics- 
based model that is derived from first principles, therefore, 
model-based approaches can provide EOL and RUL estimates 
that are much more accurate and precise than data-driven ap- 
proaches, if the models are correct. 

In this paper, we develop a general model-based prognostics 
framework using particle filters, based on preliminary work 
presented in (Daigle & Goebel, 2009, 2010b). Particle filters 
are nonlinear state observers that approximate the posterior 
state distribution as a set of discrete, weighted samples. They 
have become a popular methodology in the context of prog- 
nostics, where they are used for joint state-parameter estima- 
tion, e.g., (Saha & Goebel, 2009; Orchard, 2007; Daigle & 
Goebel, 2009, 2010b). Although suboptimal, the advantage 
of particle filters is that they can be applied to systems which 
may be nonlinear and have non-Gaussian noise terms, where 
optimal solutions are unavailable or intractable. Further, be- 
cause they are based on probability distributions, they help 
in managing the uncertainty that may arise from a number of 
sources in prognostics. 

As a case study, we consider a pneumatic valve from the 
Space Shuttle cryogenic refueling system. We construct a 
detailed physics-based model of a pneumatic valve that in- 
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eludes models of different damage mechanisms. We run a 
comprehensive set of prognostics experiments in simulation 
to demonstrate the approach and establish that prognostics 
may be performed for pneumatic valves using only discrete 
position sensors. Further, we demonstrate the approach using 
historical pneumatic valve data from the refueling system. 
The paper is organized as follows. Section 2 formally defines 
the prognostics problem and describes the computational ar- 
chitecture. Section 3 presents the modeling methodology and 
develops the model of the pneumatic valve. Section 4 dis- 
cusses the damage estimation approach using particle filters, 
and Section 5 provides the prediction algorithm. Section 6 
presents simulation results and the demonstration of the ap- 
proach on real data. Section 7 discusses related work. Sec- 
tion 8 concludes the paper. 

2 Prognostics Approach 

The problem of prognostics is predicting the EOL and/or the 
RUL of a component. In this section, we first formally define 
the problem of model-based prognostics. We then describe a 
general model-based architecture within which a prognostics 
solution may be implemented. 

2.1 Problem Formulation 

We assume the system may be described by 

x(f) = f (t, x(t), 0(t),u(t),v(t)) 

y(t) = h(t, x(f), 9{t), u (t), n (t)) 

where f £ R is the continuous time variable, x(t) £ R”* 
is the state vector, 9(t) £ R" 8 is the parameter vector, 
u(t) £ R" 1 * is the input vector, v(t) £ R™ 1 ’ is the process 
noise vector, f is the state equation, y (t) £ R nH is the output 
vector, n(t) £ R"" is the measurement noise vector, and h is 
the output equation. This representation considers a general 
nonlinear model with no restrictions on the functional forms 
of f or h. Further, the noise terms may be coupled in a nonlin- 
ear way with the states and parameters. The parameters 9(t) 
evolve in an unknown way. In practice, they are typically 
considered to be constant, but may in fact be time-varying. 
Our goal is to predict EOL at a given time point t p using the 
discrete sequence of observations up to time tp, denoted as 
yo . tp . EOL is defined as the time point at which the compo- 
nent no longer meets one of a set of functional requirements 
(e.g., a valve does not open in the required amount of time). 
These requirements may be expressed using a threshold, be- 
yond which we say the component has failed. In general, 
we may express this threshold as a function of the system 
state and parameters, Teol( x (/), 9(t)), which determines 
whether the system has failed, where Teol(x(/), 9(t)) = 1 
if a requirement is violated, and 0 otherwise. Using this func- 
tion, we can formally define EOL with 

EOL(tp) = inf{f £ R : t > t P A Teol( x (/), 0(f)) = 1}, 


and RUL with 

RUL{t P ) = EOL{t P ) - t P . 

In practice, many sources of uncertainty exist that unfold into 
the prediction. Uncertainty in the initial state, uncertainty in 
the model, process noise (i.e., v(f)), and sensor noise (i.e., 
n(f)) result in an uncertain estimate of (x(f), 0(f)). In pre- 
dicting from this uncertain state, modeling errors, process 
noise, and uncertainty in the future inputs of the system fur- 
ther add to the overall uncertainty. At best, then, we can only 
compute a probability distribution of the EOL or RUL, rather 
than a single prediction point. The goal, then, is to compute, 
at time t P , p(EOL(t p )\y 0 . tp ) and/or p(RUL(t P )\y 0: t P )- 

2.2 Prognostics Architecture 

We adopt a model-based approach, wherein we develop de- 
tailed physics-based models of components and systems that 
include descriptions of how faults evolve in time. These mod- 
els depend on unknown and possibly time- varying parame- 
ters. Therefore, our solution to the prognostics problem takes 
the perspective of joint state -parameter estimation. In discrete 
time k, we estimate x/, ; and 9k, and use these estimates to pre- 
dict EOL and RUL at desired time points. 

We employ the prognostics architecture in Fig. 1. The sys- 
tem is provided with inputs ri/, and provides measured out- 
puts yfc. Prognostics may begin at k = 0, with the dam- 
age estimation module determining estimates of the states 
and unknown parameters, represented as a probability distri- 
bution p(xfc, 9k |yo:fc). The prediction module uses the joint 
state -parameter distribution, along with hypothesized future 
inputs, to compute EOL and RUL as probability distributions 
p(EOL kp |y 0:fcp ) and p{RUL kp \y 0 :k P ) at given prediction 
times kp. In parallel, a fault detection, isolation, and identifi- 
cation (FDII) module may be used to determine which dam- 
age mechanisms are active, represented as a fault set F. The 
damage estimation module may then use this result to limit 
the space of parameters that must be estimated. Alternatively, 
prognostics may begin only when diagnostics has completed. 
In this paper, we assume prognostics is not aided by FDII. 

3 Modeling Methodology 

To implement a model-based prognostics approach, detailed 
physics-based models of the component are necessary. Such 
models must not only describe the nominal behavior of the 
component, but also the faulty behavior. Further, they must 
describe how faults, or damage, progress in time. It is with 
these models that prediction may be performed. 

To construct models in this paradigm, the first step is to de- 
velop a nominal model based on a first principles, physical 
understanding of the system. The next step is to identify the 
variables or parameters that are representative of faults, which 
we call damage variables, d(f), e.g., a friction parameter, or 
the size of a leak. Since these variables change dynamically. 
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Figure 1. Prognostics architecture 



they form part of the state vector, i.e., d(t) C x(t). The final 
step is to develop models for how these variables change in 
time, i.e., how the damage progresses or evolves, which we 
call damage progression functions. These models are typi- 
cally dependent on other states of the system, and, therefore, 
dependent on the system inputs. They augment the state equa- 
tion f . These functions are parameterized by unknown param- 
eters, which we call wear parameters, w (f) C e(t). 

In the remainder of this section, we apply this modeling 
framework to a pneumatic valve, which serves as the case 
study for this paper. 

3.1 Pneumatic Valve Modeling 

Pneumatic valves are gas-actuated valves used in many do- 
mains. A normally-closed valve with a linear cylinder actua- 
tor is illustrated in Fig. 2. The valve is opened by filling the 
chamber below the piston with gas up to the supply pressure, 
and evacuating the chamber above the piston down to atmo- 
spheric pressure. The valve is closed by filling the chamber 
above the piston, and evacuating the chamber below the pis- 
ton. The return spring ensures that when pressure is lost, the 
valve will close due to the force exerted by the return spring, 
hence it is a normally-closed valve. 

We develop a physics model of the valve based on mass and 
energy balances. The system state for the nominal model 


includes the position of the valve, x(t), the velocity of the 
valve, v(t), the mass of the gas in the volume above the pis- 
ton, rrit{t), and the mass of the gas in the volume below the 
piston, m b (t): 

x(t) = [x(t) v(t) m t (t) m b (t)] T . 

The position is defined as x = 0 when the valve is fully 
closed. The stroke length of the valve is denoted by L s ; when 
the valve is fully open its position x = L s . 

The derivatives of the states are described by 

x(t)=[v{t) a(t) f t (t) f b {t)] T , 

where a(t) is the valve acceleration, and f t {t) and f b (t) are 
the mass flows going into the top and bottom pneumatic ports, 
respectively. 

The inputs are considered to be 

u (f) = [pi(t) p r (t) U t (t ) U b (t)\ T , 

where Pi(t) and p r (t) are the fluid pressures on the left and 
right side of the plug, respectively, and u t (t) and u b (t) are the 
input pressures to the top and bottom pneumatic ports. These 
pressures will alternate between the supply pressure and at- 
mospheric pressure depending on the commanded valve posi- 
tion. 

The acceleration is defined by the combined mass of the 
piston and plug, m, and the sum of forces acting on the 
valve, which includes the forces from the pneumatic gas, 
(p b (t) — p t (t))A p , where p b (t) and p t (i) are the gas pres- 
sures on the bottom and the top of the piston, respectively, 
and A p is the surface area of the piston; the forces from the 
fluid flowing through the valve, (p r (t) — pi(t))A v , where A v 
is the area of the valve contacting the fluid; the weight of the 
moving parts of the valve, —mg, where g is the acceleration 
due to gravity; the spring force, —k(x(t) — x a ), where k is 
the spring constant and x a is the amount of spring compres- 
sion when the valve is closed; friction, —rv(t), where r is the 
coefficient of kinetic friction, and the contact forces F c (t) at 
the boundaries of the valve motion, 

! k c (—x), if x < 0, 

0, if 0 < a; < L s , 

-k c (x-L s ), if x > L s , 
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where k c is the (large) spring constant associated with the 
flexible seals. Overall, the acceleration term is defined by 


a{t) = — \(Pb(t) — p t (t))A p + {p r (t) - pi(t))A v 

TO L 

— mg — k(x(t) — x Q ) — rv(t) + F c (t ) 


The pressures pt(t) and pb(t) are calculated as: 

m t (t)R g T 
Vt 0 + A P (L S — x(t)) 
m b (t)R g T 


Pt{t ) = 

Pb(t) = 


Vb 0 + A p x(t) 


where we assume an isothermal process in which the (ideal) 
gas temperature is constant at T, R g is the gas constant for the 
pneumatic gas, and l) n and V bo are the minimum gas volumes 
for the gas chambers above and below the piston, respectively. 
The gas flows are given by: 


ft(t) = fg(pt(t),u t {t)) 

fb(t ) = fg(p b (t),U b (t)), 


where f g defines gas flow through an orifice for choked and 
non-choked flow conditions (Perry & Green, 2007). Non- 
choked flow for p, > p 2 is given by f g ,nc{Pi,P 2 ) = 


C s A sPl 




( El 

\Pi 


2 

7 



where 7 is the ratio of specific heats, Z is the gas compress- 
ibility factor, C s is the flow coefficient, and A s is the orifice 
area. Choked flow for pi > p 2 is given by 


/ S ,c(Pl,P2) 


C s A sP i 


\ 


7 

( 2 \ 

ZR.gT 



7+1 

-7—I 


Choked flow occurs when the upstream to downstream pres- 
sure ratio exceeds (-^) 7 ^ 7 The overall gas flow equa- 
tion is then given by 


fg(Pl,P2) = { 


(fg ,nc (Pl,P2) ifPl>P2 

and ^ <(1+1)^, 

fg,c(Pl ? P 2 ) ifPl>P2 

and > ( 2 + 1 ) 1 ^ 
P2 — v 2 J 

f g,nc (P2,Pl) if P2 > Pi 

«»<! “ < , 
-fg, c(j>2,Pl) rfP2>Pl 

and £>(¥)'*’. 


We select our complete measurement vector as 

y(t) = [x(t) p t (t) p b (t) f v (t ) open(t ) closed{t)] T . 


Here, f v is the fluid flow through the valve, given by 

fv(t) = '—C v A vS J'± p | Pfi — P/r|sign(p/( — Pfr), 

where C v is the (dimensionless) flow coefficient of the valve, 
p is the liquid density, and we assume a linear flow charac- 
teristic for the valve. The open(t ) and closed(t ) signals are 
from discrete sensors that output 1 if the valve is in the fully 
open or fully closed state: 

open(t) = 
closed(t ) = 

Fig. 3 shows a nominal valve cycle. The valve is commanded 
to open at 0 s. The top pneumatic port opens to atmosphere 
and the bottom port opens to the supply pressure (approxi- 
mately 5.3 MPa, or 750 psig). When the force on the un- 
derside of the piston is large enough to overcome the return 
spring, friction, and the gas force on the top of the piston, the 
valve begins to move upward as the pneumatic gas continues 
to flow into and out of the valve actuator. At about 8 s, the 
valve is completely open. The valve is commanded to close 
at 15 s. The bottom pneumatic port opens to atmosphere and 
the top port opens to the supply pressure. When the force bal- 
ance becomes negative, the valve starts to move downward, 
and completely closes at around 20 s. The valve closes faster 
than it opens due to the return spring. 


if x(t) > L s 

otherwise 

if x(t) < 0 

otherwise. 


3.2 Damage Modeling 

With the nominal model defined, we may now determine 
which damage mechanisms are relevant, and what model pa- 
rameters change as a result of the damage. From valve docu- 
mentation and historical maintenance records, we have iden- 
tified a set of faults that, although not exhaustive, contains 
the most important and most often observed faults that affect 
valve functionality and lead to EOL. The set of faults includes 
friction damage, spring damage, internal valve leaks, and ex- 
ternal valve leaks. The functional requirements on the valve 
that define Teol are that it opens and closes within given 
timing limits, and that it fully closes upon loss of actuating 
pressure. In this example, we use 15 s for both the opening 
and closing time limits. 

One damage mechanism present in valves is sliding wear. 
The equation for sliding wear takes on the following 
form (Hutchings, 1992): 

V(t) = w\F(t)v(t)\, 

where V(t) is the wear volume, w is the wear coefficient 
(which depends on material properties such as hardness), 
F(t) is the sliding force, and v(t) is the sliding velocity. 
Friction will increase linearly with sliding wear, because the 
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Figure 3. Nominal valve operation 
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contact area between the sliding bodies becomes greater as 
surface asperities wear down (Hutchings, 1992). Lubrication 
between the sliding bodies can also degrade over time. We 
therefore characterize friction damage as change in the fric- 
tion coefficient, and model the damage progression in a form 
similar to sliding wear: 

f{t) = w r \F f (t)v(t)\ 

where w r is the wear coefficient, and Ff(t) is the friction 
force defined in the previous subsection. The friction coeffi- 
cient only grows when the valve is moving, as only then can 
sliding wear occur. Therefore, the friction coefficient evolves 
in a step-wise fashion, with damage only growing during the 
opening and closing motions. As the friction coefficient in- 
creases, the friction force increases, further increasing the rate 
at which the friction coefficient grows. This results in a dam- 
age progression that is similar to an exponential when viewed 
at large time scales. 

Fig. 4 shows the effect of an increase in friction on the valve 
cycle. We define r + as the maximum value of the friction pa- 
rameter at which the valve still behaves within the timing lim- 
its. Above this value, the friction force becomes large enough 
that the valve cannot open within the 15 s limit, as shown in 
Fig. 4. So, T E OL(x(t), 9(t)) = 1 if r(t) > r + . 

We assume a similar equation form for spring damage: 
k(t) = -w k \F s (t)v(t)\, 


Figure 4. Valve operation with damage 

force. The more the spring is used, the weaker it becomes, 
characterized by the change in the spring constant. Like with 
friction damage, the spring constant decreases only when the 
valve moves. As the spring becomes damaged, the spring 
force will decrease, and so the rate at which spring damage 
occurs will also decrease. 

Fig. 4 shows the effect of a decrease in the spring constant 
on the valve cycle. In normal operation, without the spring 
tending the valve to close, the valve will open faster and close 
slower. The spring can actually fail completely, and the valve 
would still open and close in time, however, the spring must 
be strong enough to close the valve against system pressure 
when the actuating pressure is lost. So, it is the loss of this 
function leads to Teol for spring damage, and we define k~ 
as the smallest value of k at which the valve will still fully 
close upon loss of supply pressure. So, T E OL(x-(t), 6(t)) = 1 
also if k(t) < k~ . 

An internal leak in the valve can appear at the seal surround- 
ing the piston as a result of sliding wear. The pneumatic gas 
is then able to flow between the volumes above and below the 
piston, decreasing the response time of the valve. We parame- 
terize this leak by its equivalent orifice area, Ai(t), described 
by 

Ai(t) = Wi\F f (t)v(t)\, 


where Wk is the spring wear coefficient and F s (t) is the spring where w L is the wear coefficient. The mass flow at the leak. 
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is computed using the gas flow equation described ear- 
lier: 


fi{t) = fg(Pt(t),Pb{t)), 

where positive flow denotes flow from the top volume to the 
bottom volume. The term is subtracted from the f t (t) equa- 
tion and added to the f b (t) equation. 

As sliding wear occurs, the leak size increases. As with the 
friction and spring damage, the internal leak only grows when 
the valve is moving. Over large time scales, the internal leak 
appears to grow linearly with the valve cycles, since the fric- 
tion force does not change much as the leak size grows. 

The presence of an internal leak makes it more difficult to ac- 
tuate the valve, because it causes gas to flow into the lower 
pressure volume that is being evacuated and out of the higher 
pressure volume that is being filled. We define A+ as the 
maximum internal leak area at which the valve opens and 
closes within the functional limits. So, Tbol( x (^)i 0(f)) = 1 
also if A t (t) > A - . Fig. 4 shows the effect of an internal leak 
on the valve cycle. 

External leaks may also form, most likely at the actuator con- 
nections to the pneumatic gas supply, due to corrosion and 
other environmental factors. Without knowledge of how the 
leak size progresses, we assume the growth of the area of the 
leak holes, A e , t (t) and A e , b (t), is linear: 

A e , t (f) = W e ,t 

Ae,b(t') — W e ,bt 

where the t and b subscripts denote a leak at the top and bot- 
tom pneumatic ports, respectively, and w e ,t and w e j, are the 
wear coefficients. The leak flows are given by 

/e,t(f) = fg(Pt{t),Patm(t)) 
fe,b(t) = fg(Pb(t),Patm(t)), 

where p a tm. is atmospheric pressure. The / e>t (f) term is sub- 
tracted from the f t (t) equation, and the f e ,b(t) term is sub- 
tracted from the fb(t) equation. 

Note that this damage progression is independent of the valve 
inputs. Fig. 4 shows the effects of top and bottom external 
leaks on the valve cycle. The effect of the formation of a 
leak at the top pneumatic port is that it becomes easier to 
open the valve but more difficult to close it. Conversely, the 
effect of a leak at the bottom pneumatic port is that it be- 
comes more difficult to open but easier to close the valve. 
We define the maximum leak hole areas at which the valve 
still opens and closes within the functional limits as A+ t and 
A+ b . So, T E oL{x(t), 0(t)) = 1 also if A e , t (t) > or 
A e ,b(t) > A, h . Alternatively, maximum allowable leakage 
rates may define EOL. 

So, the damage variables are given by 

d(t) = \r(t) k(t) Ai(t) A e , t (t) A e ^ b (t)] T , 


and the complete state vector of the extended valve model 
becomes 


x(f) 


x(t)v(t) 
m t (t) 
m b {t ) 
r(t) 
k(t) 
Ai{t ) 

A e , b (t ) 


The wear parameters form the unknown parameter vector, i.e.. 


w(i) = 6(t) 


w r (t) 

w k {t) 

m(t) 

W e , bit ) 


4 Damage Estimation 


In the model-based paradigm, damage estimation reduces 
to joint state-parameter estimation, i.e., computation of 
p(xfc, 0fc|yo:fc). The typical approach to state estimation is to 
use a state observer, which helps to compensate for the effects 
of process and sensor noise on the state estimate. Joint state- 
parameter estimation is traditionally performed using a state 
observer by augmenting the state vector with the parameter 
vector. In this way, both states and parameters are simultane- 
ously estimated. 

For linear systems with additive Gaussian noise terms, the 
Kalman filter is appropriate. For nonlinear systems with ad- 
ditive Gaussian noise terms, the extended Kalman filter or 
unscented Kalman filter (Julier & Uhlmann, 1997) may be 
used. However, for nonlinear systems with non-Gaussian 
noise terms, particle filters are best suited, and offer approx- 
imate (suboptimal) solutions to the state estimation problem 
for such systems where optimal solutions are unavailable or 
intractable (Arulampalam, Masked, Gordon, & Clapp, 2002; 
Cappe, Godsill, & Moulines, 2007). In particle filters, the 
state distribution is approximated by a set of discrete weighted 
samples, called particles. As the number of particles is in- 
creased, accuracy increases and the optimal solution is ap- 
proached. In addition, particle filters are straightforward to 
implement, and computational complexity may be controlled 
by increasing or decreasing the number of particles, with re- 
spect to the desired estimation performance. Further, the dis- 
crete position sensors cannot be directly handled by other fil- 
tering algorithms. For these reasons, we select particle filters 
for our model-based prognostics framework. 

With particle filters, the particle approximation to the state 
distribution is given by 

{(x*,0‘ fc ), «,*}£*, 

where N denotes the number of particles, and for particle i, 
x'j, denotes the state estimates, 0' k denotes the parameter es- 
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Algorithm 1 SIR Filter 

Inputs: {(xjUi.ei.!), Wfc-l}£l: Ufc-i:fc,yfc 

Outputs: 

for i = 1 to N do 

ei~p(e k |ei_i) 

Xfc ~ p{Xk\Xk r l,0k_l, Ufe-i) 
w\ <-p(yfc|xL,0L.Ufe) 

end for 

N 

W 4 ~ wl 

i=l 

for i = 1 to N do 

wi «- wi/w 

end for 

{(xfc,0fc).«>fc}£Li «- Resample({(x^,0 I fc ),Wfc}iI 1 ) 


timates, and tujj. denotes the weight. The posterior density is 
approximated by 

N 

J?(x fc ,0fc |y 0 :fe) « ^wlS^ ^idxkdGk), 

2=1 

where <J/ X * d i\(dxkdOk ) denotes the Dirac delta function lo- 
cated at (x|, 0 l fc ). 

We employ the sampling importance resampling (SIR) parti- 
cle filter, and implement the resampling step using systematic 
resampling (Kitagawa, 1996). The pseudocode for a single 
step of the SIR filter is shown as Algorithm 1. Each particle i 
is propagated forward to time k by first sampling new parame- 
ter values. Here, the parameters Of, evolve by some unknown 
random process that is independent of the state x/ ;: . To per- 
form parameter estimation within a particle filter framework, 
however, we need to assign some type of evolution to the pa- 
rameters. The typical solution is to use a random walk, i.e., 
for parameter 0, 9 *. = 0fc_i + £/._ i, where (,k-i is typically 
Gaussian noise. After sampling parameter values from the 
selected distribution, new states are sampled by applying the 
state equation f to (x^, 0),) with process noise v(i) sampled 
from its assumed distribution. The particle weight is assigned 
using y j~. Specifically, the output equation h is applied to 
(Xfc+i) 0fc+i) and the likelihood of the corresponding output 
is computed using the assumed probability density function of 
the sensor noise. The weights are then normalized, followed 
by the resampling step. 

Pseudocode for the resampling step is given as Algorithm 2. 
First, the cumulative distribution function (CDF) of the 
weights is computed as c. A starting point is drawn from 
the uniform distribution. The algorithm then moves along the 
CDF and a new particle is resampled from the old particle set. 
The weights of the resampled particles are all assigned to be 
the same. The resampling step is necessary to avoid degener- 
acy, in which all but a few particles have negligible weight. 
If this occurs, then the state distribution is very poorly rep- 
resented, and computation is wasted on particles that do not 
contribute to the approximation. Resampling based on the 


Algorithm 2 Systematic Resampling 

Inputs: {(xJ.,0J.) ? u4} f =1 
Outputs: { (x J k ,0{), w 3 k }f =1 

C 1 4 — 0 

for i = 2 to N do 

C i 4 C;_i + Wf~ 

end for 

*4—1 

Ml ~ U(0, 1/N) 

for j = 1 to A’ do 

U 4 - Mi + (j + 1 )/N 

while m > c i do 

i 4 — i + 1 

end while 

x l x fe 

K <- ei 

wi = l/N 

end for 


CDF avoids this problem by multiplying particles with high 
weight and dropping particles with low weight. Here, we re- 
sample at each time step, but the frequency of resampling may 
be reduced by computing the effective number of particles, 
and only resampling when this number falls below a given 
threshold (Arulampalam et al., 2002). 

During the sampling step, particles are generated with param- 
eter values that will be different from the initial guesses for 
the unknown parameters. The particles with parameter val- 
ues closest to the true values should match the outputs bet- 
ter, and therefore be assigned higher weight. Resampling 
will cause more particles to be generated around these bet- 
ter values. As this process is repeated over time, the parti- 
cle filter converges to the true values. The selected variance 
of the random walk noise must be large enough so as to al- 
low convergence in a reasonable amount of time, but small 
enough such that when convergence is reached, the parame- 
ter can be tracked smoothly. Fig. 5 illustrates the effect of 
the variance value on estimation performance shown in an 
accuracy-precision plot. The horizontal axis provides esti- 
mation accuracy as computed using percent root mean square 
error (PRMSE) of the unknown wear parameter, w f .j n and 
the vertical axis provides the variability of the wear param- 
eter distribution as computed using relative median absolute 
deviation (RMAD) (performance metrics are defined in the 
Appendix). The optimal point is located at the origin. As 
the random walk variance of w e y, is decreased, both PRMSE 
and RMAD become smaller, improving in both performance 
dimensions. As the random walk variance is decreased be- 
yond 5 x 1CP 20 , however, the wear parameter estimate does 
not always converge and cannot be tracked. Since the pa- 
rameter values are unknown, an optimal value cannot be de- 
termined a priori. So, selecting an appropriate value can be 
difficult, but knowledge of the correct order of magnitude of 
the parameter is helpful. Additionally, correction loop meth- 
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Algorithm 3 EOL Prediction 


Inputs: {(xl p ,6 t k ),wl I 
Outputs: {EOU kp ,wl 
for i = 1 to A do 
k t — k p 


\ N 

L=i 


01 


01 


while Teo z, (xj. , 0J. ) = 0 do 

Predict u k 


0 


fc+i 
c fc + l 


‘P(0k+ i\0i: 
|xl 


' p(x k+1 \x.l,0i,u k ) 


k <— k + 1 


01 


0 


'-fc+i 

i 

fc+1 


end while 

EOL{ p i k 

end for 


Figure 5. Estimation performance for a bottom external leak, 
with N = 500, M = {x, f,Pt,Pb }• and varying random walk 
variance 

ods can be used to tune this value online as a function of per- 
formance (Orchard, Kacprzynski, Goebel, Saha, & Vachtse- 
vanos, 2008; Saha & Goebel, 2011; Daigle & Goebel, 201 1). 
If the unknown parameters may be assumed to be constant, 
then other approaches can be employed to improve estimates 
and offset the increase in covariance contributed by the ran- 
dom walk (Liu & West, 2001; Clapp & Godsill, 1999). 

Note that the discrete position sensors ( open and closed ), in 
reality, have no noise, but a certain amount of sensor noise 
must be assumed for sensors within the particle filter frame- 
work. Therefore, within the algorithm, we assume some sen- 
sor noise is present for these sensors. 

5 Prediction 

Prediction is initiated at a given time kp. Using the cur- 
rent state estimate, p(x k ,O k \y 0 . k ) the goal is to compute 
p(EOL kp |y 0 :fe P ) and p(RUL kp |y 0;fcp ). The particle filter 
computes 

N 

p(x kp , e kp |y 0:fep ) « ^2 w l kp <5 (x » ^ ^ ) (dx kp dd kp ). 

i= 1 

We can approximate a prediction distribution n steps forward 
as (Doucet, Godsill, & Andrieu, 2000) 

Pfa-kp+m 0kp+n\yO:kp') ~ 

N 

^ ] W kp\yt i kp + n ,O i kp+rl )(d x kp+nd0 k p-\-n)- 
i—1 

So, for a particle i propagated a steps forward without new 
data, we can simply take its weight as w kp . Similarly, we can 


approximate the EOL as 

N 

p(EOL kp \y 0:kp ) « Y,<P*EOL ip (dEOL kp ). 

i—l 

To compute EOL, then, we propagate each particle forward 
to its own EOL and use that particle’s weight at kp for the 
weight of its EOL prediction. 

The pseudocode for the prediction procedure is given as Al- 
gorithm 3 (Daigle & Goebel, 2010b). Each particle i is propa- 
gated forward until Teol(^\, 0 k ) evaluates to 1; at this point 
EOL has been reached for this particle. The complexity of 
the algorithm is variable, as each particle must be simulated 
forward until its individual EOL, which is different for each 
particle. 

Prediction requires hypothesizing future inputs of the system, 
U/g, because damage progression is rarely independent of the 
system inputs. In general, these inputs are unknown and must 
be chosen in order to satisfy both operational constraints and 
the type of prediction that is required. For the valve, we as- 
sume input-independence for the external leaks, as described 
by their damage models, but not for the remaining faults. For 
the valve, the problem of hypothesizing inputs is simplified. 
Each valve cycle corresponds to the same set of inputs, except 
for the fluid pressures p k and pp. However, these pressures 
can safely be assumed to be constant in our application do- 
main, and, further, they have an almost negligible effect on 
valve behavior because the forces they produce are very small 
compared to the other forces acting on the valve. So, the fu- 
ture inputs for each valve cycle are deterministic. We can 
simply provide repeated valve cycles as input, and the pre- 
diction step will determine after how many valve cycles EOL 
will be reached for each particle. 

Fig. 6 shows results from the prediction of friction damage for 
N = 100. Initially, the particles have a very tight distribution 
of friction parameter values, but the distribution of the wear 
parameter, w r , is relatively large. As a result, the individual 
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Parameter 

Value 

9 

9.8 m/s 2 

Patm 

1.01 x 10 5 Pa 

P supply 

5.27 x 10 6 Pa 

m 

50 kg 

r 

6.00 x 10 3 Ns/m 

k 

4.80 x 10 4 N/m 

X 0 

2.54 x 10 _1 m 

L s 

3.81 x 10~ 2 m 

A p 

8.10 x 10~ 3 m 2 

Vto 

8.11 x 10“ 4 m 3 

v b0 

8.11 x 10~ 4 m 3 

P 

7.10 x 10 2 kg/m 3 

A v 

5.07 x 10" 2 m 2 

C v 

4.36 x 10 _1 

Rg 

2.96 x 10 2 J/K/kg 

T 

293 K 

7 

1.4 

Z 

1 

A s 

1.00 x 10~ 5 m 2 

C s 

0.62 


Table 1 . Nominal valve parameters 


Figure 6. Prediction of friction damage 

trajectories for the different wear parameter values are eas- 
ily distinguishable as EOL is approached. The different EOL 
values along with particle weights form an EOL distribution 
approximated by the probability mass function shown in the 
figure. 

6 Results 

We apply the prognostics framework to a pneumatic valve of 
the Space Shuttle cryogenic refueling system. This system 
transfers cryogenic propellant (liquid hydrogen) from a stor- 
age tank to the vehicle tank through a network of pipes and 
valves. We focus on one of the transfer line valves in this sys- 
tem. In this section, we first develop and validate the model 
for the selected valve. Then, since only the discrete open 
and closed sensors are available for the valve, we establish, 
through a set of simulation-based experiments, that prognos- 
tics can still be performed using only those sensors. We then 
apply the prognostics approach to historical valve degradation 
data. Since proactive maintenance was performed early, no 
complete run-to-failure data exist. We augmented the miss- 
ing portion with simulated data. 

6.1 Model Validation 

For the selected pneumatic valve, we identified model param- 
eters using component specifications and estimated the re- 
maining parameters (to, r, k, x Q , and C s ) from nominal data. 
The model parameter values are shown in Table 1 . Note that 
the pneumatic gas is nitrogen. 


Nominally, this valve opens in a little over 1 s and closes 
within 0.2 s. The valve has an opening time requirement 
of 5 s and a closing time requirement of 4 s, and these de- 
fine functional thresholds for EOL. The valve is considered to 
have failed when these timing limits are exceeded. Further, 
the valve must fully close when actuating pressure is lost. 

Due to the high safety margin of the fueling system, com- 
ponents must always meet tight operational constraints, and 
components are never intentionally run to failure. We focus 
on a particular time frame of historical data, consisting of 7 
valve cycles, where a clear trend in performance degradation 
was observed before a maintenance action took place. We an- 
alyze this particular data set to determine the dominant dam- 
age mode and validate the corresponding damage model. 

The obtained valve data includes only the open and closed 
sensors. Along with valve open/close commands, we may 
extract the following timing information, shown in Fig. 7: 

1. t open ,i, defined as the difference between the time when 
the valve is commanded to open and it starts to move 
open, 

2. topen, 2 , defined as the difference between the time when 
the valve starts to move open and it fully opens, 

3- t open = i 0 pera,i + topen, 2 - which is the total time for the 
valve to open, 

4. trioseA , defined as the difference between the time when 
the valve is commanded to close and it starts to move 
closed, 

5- fciose, 2 , defined as the difference between the time when 
the valve starts to move closed and it fully closes. 
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Figure 7. Pneumatic valve timing diagram. A command of 1 
opens the valve, and a command of 0 closes the valve. 
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Figure 8. Pneumatic valve model comparison for growth of 
bottom external leak 


6. tclose = tclose, 1 + t c i a se ,2, which is the total time for the 
valve to close. 

We plot these values for the degrading valve in Fig. 8. We 
note trends in each of the timing variables. Most importantly, 
the opening time increases, while the closing time decreases. 
From Fig. 4 (Section 3), we know that only two fault modes 
may cause this behavior, namely, the internal leak, and the 
bottom external leak. We may distinguish between these two 
damage modes by looking at t npfm ,\ ■ For the internal leak, 
a slight decrease in this value is expected as the damage pro- 
gresses, but for the bottom external leak, a significant increase 
in this value is expected. According to the data, this value 
increases, therefore, the bottom external leak is the only con- 
sistent damage mode, and must be dominant. Fig. 8 shows 
simulated timing values plotted alongside the measured tim- 
ing values, with a wear rate of w e ^ = 3.0 x 10“ ' m 2 /cycle, 
which best matched the data. The simulated data captures all 
the trends of the real data, and the simulated timing values are 
all fairly close to the real values at each cycle. 


6.2 Measurement Analysis 

Since only the discrete open and closed sensors are avail- 
able, we must first establish that prognostics can still be per- 
formed satisfactorily in this case. To investigate this issue, 
we performed a number of simulation experiments for each 
fault under different measurement sets. In each case, we used 
N = 500 particles, and tuned the random walk variances 
assuming that the order of magnitude of the wear parame- 
ters was known. Wear parameter values were selected so that 
EOL occurred near 100 cycles. For each fault and measure- 
ment set, 5 experiments were performed. Table 2 presents 
averaged results over these experiments. The performance 
metrics are defined in the Appendix. Estimation is evalu- 
ated based on the estimate of the unknown wear parameter, 
as it is this estimate that most greatly influences subsequent 
prognostics performance, since the future valve usage is well- 
defined. Estimation accuracy is calculated using the percent- 
age root mean square error (PRMSE), where we ignore the 
first 10 cycles, which are associated with convergence. We 
calculate the spread using the relative median absolute devi- 
ation (RMAD), which is a robust measure of statistical dis- 
persion. We compute also convergence of the wear parameter 
estimate, denoted as C w , over the first 3 cycles. For a pre- 
diction point kp, we compute measures of accuracy and pre- 
cision. For accuracy, we use the relative accuracy (RA) met- 
ric (Saxena, Celaya, Saha, Saha, & Goebel, 2010). We use 
RA to denote the RA averaged over all prediction points. We 
calculate prediction spread using RMAD, which we denote 
as RMAD/jjtx. We use RMAD/.f-/ to denote R M A D RU L 
averaged over all prediction points. 

From Table 2, it is clear that prognostics can be performed 
successfully with all measurement sets, including using only 
the discrete sensors open and closed. Because of the loss of 
information in going from the continuous sensors to the dis- 
crete sensors, the state distribution has a much wider variance, 
usually increasing by around 50%. In turn, the wider variance 
causes the estimate to smooth out, and, as a result, PRMSE 
is lowered, and this corresponds to increases in RA, usually 
around 0.5%, however this is at the cost of a less confident 
prediction, as shown by RMAD^j/^. Further, convergence of 
the wear parameter estimate is slower using only the discrete 
sensors. 

The results show that, if it is feasible to add the continuous 
sensors, prognostics performance can be improved signifi- 
cantly. The increase in information provided by these sen- 
sors leads to more confident predictions. Different sensors 
are more useful for different faults, so additional sensors may 
be selected based on which faults are more likely to appear, 
or more critical. For friction and spring damage, the posi- 
tion sensor is more useful than the pressure sensors, as shown 
by the decreased spread when M = {a;} compared to when 
M = {pt,Pb} for these faults. In contrast, the pressure sen- 
sors are more valuable for the leak faults. 
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Fault 

M 

PRMSE m 

RMAD W 

c w 

RAiiyi. 

RM AD /((//, 

Friction 

{x,f,pt,p b } 

4.00 

8.07 

30.72 

96.59 

8.64 


{*} 

3.92 

9.38 

28.79 

97.24 

8.62 


{Pt,Pb} 

3.58 

10.11 

34.90 

96.27 

10.17 


{open, closed} 

3.52 

13.53 

44.68 

97.38 

11.68 

Spring Rate 

{x,f,pt,p b } 

5.58 

10.87 

32.70 

94.91 

11.04 


{*} 

5.84 

12.36 

36.07 

94.80 

12.43 


{Pt,Pb} 

5.47 

13.14 

31.74 

94.51 

13.09 


{open, closed} 

5.10 

16.58 

43.50 

94.58 

16.40 

Internal Leak 

{x,f,pt,p b } 

6.00 

10.33 

25.88 

95.26 

10.82 


W 

5.84 

13.81 

21.40 

95.29 

14.37 


{Pt,Pb} 

5.39 

10.52 

27.59 

94.74 

10.69 


{open, closed} 

4.51 

17.42 

33.49 

95.94 

17.06 

Top External Leak 

{x,f,pt,p b } 

5.33 

9.27 

21.54 

95.73 

9.32 


{x} 

6.78 

10.99 

17.30 

94.72 

11.06 


{Pt,Pb} 

4.00 

10.04 

20.83 

96.87 

10.12 


{open, closed} 

5.14 

14.11 

18.82 

95.94 

14.24 

Bottom External Leak 

{x,f,p t ,Pb} 

5.52 

11.21 

25.29 

95.57 

11.13 


{z} 

6.02 

14.49 

25.15 

95.04 

14.22 


{Pt,Pb} 

5.13 

11.36 

24.45 

95.64 

11.27 


{open, closed} 

5.07 

18.59 

34.12 

95.69 

18.14 


Table 2. Prognostics performance for different measurement sets 


The valve considered here only sees discrete operation, where 
the valve is either fully closed or fully opened. In this case, 
the discrete open and closed sensors provide enough informa- 
tion for prognostics. For situations where the valve’s position 
is controlled in a continuous manner, a continuous position 
sensor is necessary for prognostics. For such valves, this sen- 
sor is usually already available for feedback to the position 
controller. 

6.3 Demonstration of the Approach 

Given that prognostics may still be performed using only the 
open and dosed sensors, we now apply the prognostic algo- 
rithms to the historical data. Because the actual data con- 
sisted of only 7 valve cycles before maintenance was per- 
formed, we extend the data with simulated valve data, based 
on the identified bottom external leak damage mode with 
We.b = 3.0 x 10“ ' m 2 /cycle, up to EOL at 38 cycles, at which 
the valve fails to open within the 5 second limit. Fig. 9 shows 
the valve timing data used for the demonstration. It is impor- 
tant to note that, even though the trend in valve opening time 
appears linear over the first 10 to 15 cycles, the overall trend 
is clearly nonlinear. This behavior is captured by the model, 
and serves as a justification of a model-based approach as op- 
posed to simple trending strategies. Projecting a linear trend 
out from the 10th cycle would produce a severe overestimate 
of the true RUL. 

To demonstrate a general solution, we allowed the particle 




Figure 9. Real and simulated valve timing 
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Figure 10. Estimation of pneumatic valve damage modes 


filter to jointly estimate all the damage modes of the valve. 
Because of the large state space, N = 2000 was used. The 
algorithm should converge to correct wear rates for each dam- 
age mode, and, in this case, identify the bottom external leak 
as the dominant damage mechanism. Estimation results are 
shown in Fig. 10. The algorithm correctly identifies the bot- 
tom external leak as the dominant damage mode, since it is 
the only damage mode for which the estimated damage vari- 
able is significantly increasing. During the real data portion, 
the algorithm is still converging, and by the time the real data 
portion has ended, the algorithm estimates the wear parameter 
within 30% of the true value, and the mean RUL prediction 
is within 25% of the true RUL. The estimated value of the 
RUL at that time point is less than the true value, resulting in 
a conservative estimate. By 10 cycles, the mean RUL predic- 
tion has improved to within 7% of the true RUL. 

The prognostics performance is summarized in Table 3. All 
times are given in cycles. RUL* denotes the true RUL value. 


We show the predictions starting at 5 cycles. At this point, 
the particle filter is still converging, but the RUL is predicted 
with an accuracy of 80%. Afterwards, RA ranges from 87 
to 97%. The RMAD of the prediction varies from around 5 
to 7%. So, RUL is predicted with reasonable accuracy, and 
the predictions are fairly confident. It is also useful to inspect 
the RUL prediction at a selected level of confidence. Fig. 11 
shows the predictions at the 99% confidence level, i.e., the 
value at which 99% of the prediction distribution is greater 
than or equal to that value. The predictions at this confidence 
level always lie below the true RUL, therefore, they serve as 
conservative estimates upon which risk-averse decisions can 
be made. 

Although encouraging, the results presented here represent 
only a limited validation of the overall approach, since the 
historical data covers only a single fault mode of the five con- 
sidered, and extends only to 7 valve cycles. A more complete 
validation requires a more flexible testbed, such as one that 
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tp 

RUL* 

RUL 

RA L 

R MAI) /,•(//. 

5.0 

33.00 

26.60 

80.60 

7.28 

7.5 

30.50 

27.53 

90.27 

5.49 

10.0 

28.00 

26.04 

92.99 

7.65 

12.5 

25.50 

24.42 

95.77 

5.96 

15.0 

23.00 

21.75 

94.57 

7.61 

17.5 

20.50 

21.84 

93.46 

6.04 

20.0 

18.00 

16.55 

91.93 

6.52 

22.5 

15.50 

13.53 

87.29 

6.63 

25.0 

13.00 

11.64 

89.52 

7.87 

27.5 

10.50 

10.15 

96.71 

6.90 

30.0 

8.00 

6.78 

84.74 

5.14 

32.5 

5.50 

5.34 

97.11 

4.81 

35.0 

3.00 

2.93 

97.69 

7.16 


Table 3. Pneumatic valve prognostics performance 
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Figure 11. RUL predictions at the 99% confidence level 


would allow the injection of faults and/or valves to be run to 
failure. 


Model-based prognostic techniques have been applied to 
other systems, and particle filtering techniques have been par- 
ticularly popular. In (Saha & Goebel, 2009), a particle filter- 
ing approach is applied to end of discharge and end of life 
prediction for lithium-ion batteries. In (Orchard et al., 2008), 
a particle filter-based diagnosis and prognosis framework is 
applied to crack growth in aircraft components. They also im- 
plement outer correction loops to control the variance of the 
random walk of the unknown parameters based on prediction 
error. This type of approach is preferred to the fixed vari- 
ance method employed in this paper, however, the approach 
presented in (Orchard et al., 2008) is only applicable when 
a feature can be derived that is a direct function of a single 
fault dimension, which is not the case with valves. In (Abbas, 
Ferri, Orchard, & Vachtsevanos, 2007), the authors apply a 
particle filter-based prognosis method to prediction of battery 
grid corrosion. A physics-based methodology to prediction of 
aircraft engine bearing spalls is developed in (Bolander, Qiu, 
Eklund, Hindle, & Rosenfeld, 2010). Particle filters are used 
for estimation of spall size, using diagnostic information from 
vibration and oil debris sensors in the update step. Prediction 
is performed using anticipated future load and speed, and the 
approach was evaluated in a physical testbed. Although the 
size of our state-parameter vector is much larger than those 
considered in these previous works, our particle filter-based 
approach does not differ from these previous approaches in 
any fundamental way. 

Other filtering techniques may also be used. In (Chelidze, 
2002), the system is modeled in a hierarchical fashion, with a 
fast-time observable subsystem coupled to a slow-time dam- 
age subsystem, where damage is tracked using a Kalman fil- 
ter. The approach is applied to tracking open circuit battery 
voltages in a vibrating beam testbed using only beam dis- 
placement measurements. In a related approach, a model- 
based prognosis methodology is developed in (Luo et al., 
2008) using an interacting multiple model filter for state- 
parameter estimation and prediction. The approach is applied 
to prognosis of fatigue cracks in a simulation of an automotive 
suspension system. 


7 Related Work 

Very little work has been done in valve prognostics. A health 
monitoring application for a pneumatic valve in a pressure 
control system is considered in (Gomes, Ferreira, Cabral, 
Glavao, & Yoneyama, 2010). Friction damage, spring ag- 
ing, and internal leaks are also identified as common failure 
modes in their application. The distribution of pressure sig- 
nals is monitored using a probability integral transform, and 
a measure of dissimilarity from a baseline distribution serves 
as a health index. The approach is purely data-driven, rely- 
ing on a nominal distribution of pressure signals to serve as 
a baseline, and comparing future pressure distributions to the 
baseline. No identification of the form of damage present is 
performed, and prediction is not performed either. 


8 Conclusions 

In this paper, we developed a model-based prognostics ap- 
proach using particle filters for joint state-parameter estima- 
tion. The estimated damage state of the valve is propagated 
forward in time to predict EOL and form EOL and RUL dis- 
tributions. We applied the approach to a pneumatic valve, de- 
veloping a detailed physics-based model that included dam- 
age progression models. Through simulation experiments, we 
established that prognostics could be performed successfully 
using only the discrete position sensors of the valve, and ap- 
plied the approach to real valve data under this constraint. The 
results here demonstrated the effectiveness of a model-based 
approach, and gave insight into the selection of sensors for 
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valve prognostics. 

Here, we advocate a model-based approach, where perfor- 
mance will depend greatly on the accuracy of the model 
provided. Clearly, developing such models is a key obsta- 
cle. The model developed here is suitable for a normally- 
closed pneumatic valve with a piston-based actuation mech- 
anism. With small modifications, the model may also be 
used for normally-open valves or valves that use a diaphragm 
mechanism. With other valve types, the actuation mecha- 
nism changes, so a new model must be developed. How- 
ever, there will be several commonalities between the dif- 
ferent models. For instance, in past work, we have applied 
similar modeling techniques to a solenoid valve (Daigle & 
Goebel, 2010a). Here, damage modes for the spring and fric- 
tion are also present, and the same damage progression mod- 
els are used. The damage models can also be translated into 
a forms for rotational motion, as in (Daigle & Goebel, 201 1), 
where a modified friction damage progression is used in a 
centrifugal pump model. Of course, one must also consider 
the cost of developing such detailed models against the cor- 
responding prognostics performance achieved, and this is the 
subject of recent work (Daigle et al., 2011). 

The valves considered here are always operated in the same 
way; they fully open, then fully close. Therefore, future valve 
inputs can be easily hypothesized, and in this paper we as- 
sumed a single known future input trajectory to provide EOL 
and RUL in the number of valve cycles. In general, the future 
inputs may not be well-known, and different sets or distribu- 
tions of possible inputs must be considered. In future work, 
we will extend the approach to these cases. Further, although 
the algorithm assumes all damage modes may be progressing 
simultaneously, the results here cover cases where only a sin- 
gle damage mode is progressing. Recent work considers the 
case when multiple damage modes are active for a centrifugal 
pump simulation (Daigle & Goebel, 2011), and investigat- 
ing this for pneumatic valves, especially with limited sensing, 
where fault masking becomes an issue, is considered as fu- 
ture work. Applying the approach to additional testbeds is 
also part of future work. 


Appendix 


Performance Metrics 


Estimation performance is evaluated based on the estimate of 
the unknown wear parameter. Estimation accuracy is calcu- 
lated using the percentage root mean square error (PRMSE), 
which expresses relative estimation accuracy as a percentage: 


PRMSE,,, = 100 


\ 


Mean/. 

f ™k-w * k \ 2 


[v K ) \ 


where u>k denotes the estimated wear parameter value at time 
k, w k denotes the true wear parameter value at k, and Mean/ 
denotes the mean over all values of k. In computing PRMSE, 


we typically ignore the initial portion of data associated with 
convergence. 

We calculate the spread using the relative median absolute 
deviation (RMAD), which expresses the spread relative to the 
median as a percentage: 


RMAD (A) = 100 


Median^ (|Aj — Medianj(X,-)|) 
Median j (A,- ) 


where X is a data set and A; is an element of that set. For es- 
timation spread, for time k we use the distribution of wear pa- 
rameter values given by the particle set at k as the data set. We 
summarize RMAD over an experiment by averaging RMAD 
over all k. We denote the average RMAD over multiple k 
using R.V1AD,, : 


R.V1AD,, = Me a n /, (R M A D , /„ ) , 


where RMAD,,,/. denotes the RMAD of the wear parameter 
at time k. 

The final estimation metric is convergence of the wear pa- 
rameter estimate, denoted as C w . We use the definition of the 
convergence metric described in (Saxena et al., 2010), where 
the convergence of a curve is expressed as the distance from 
the origin to the centroid under the curve (where a shorter dis- 
tance is better). We use the absolute error of the hidden pa- 
rameter estimate as the curve. We compute convergence only 
over the initial convergence period, so that errors after con- 
vergence do not contribute. The same time window must be 
used consistently over all experiments to properly compare. 
For a prediction point kp, we compute measures of accu- 
racy and precision. For accuracy, we use the relative accuracy 
(RA) metric (Saxena et al., 2010), computed as a percentage: 


RAi „ = 100 1 - 


\RUL%p — Median; (.R(7Z4, p ) | 
RULt 


where RUL kp is the true RUL at time kp. Here, we have 
chosen the median as a robust measure of central tendency 
of the distribution. We use RA to denote the averaged RA 
over all prediction points. We calculate prediction spread us- 
ing RMAD, which we denote as RMAD/,/,/. We average 
RMAD/;/-/ over all prediction points to obtain a single num- 
ber representing prediction spread over a single experiment, 
denoted using RMAD rul- 
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Nomenclature 


EOL 

end of life 

RUL 

remaining useful life 

t 

time (continuous) 

k 

time (discrete) 

tp or kp 

time of prediction 

X 

state vector 

e 

parameter vector 

V 

process noise vector 

y 

output vector 

n 

measurement noise vector 

Teol 

EOL threshold 

PRMSE 

percent root mean square error 

RMAD 

relative median absolute deviation 

RA 

relative accuracy 
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